Subroutine make_np_channel
use type
use var
use parallel
implicit none
integer :: i
integer, dimension(:,:,:), allocatable :: temp_po 

allocate (temp_po(3,int(4*radius**3),number_NP))
NPnodes = 0
NPbnodes = 0
temp_po = 0
NP = .false.
boundary_NP = .false.

call grid_NP()
call find_normal1()
if (restart==0) call anchoring_q()
call find_boundary_neighbors1()

deallocate(temp_po)
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine grid_NP()
use var
implicit none
integer(i4b) :: i,j,k, ii , label
integer(i4b) :: cont_b = 0 , cont_np = 0, R
real(dp) :: checkx, checky, checkz, check1, check2, check3
real(dp) :: Lx, Ly, Lz
real, dimension(3) :: xs

R = radius

Do ii=1,number_NP
cont_b = 0
cont_np = 0
xs = NP_position(:,ii)

if (channel) then

do k=1,dime(3)
   do j=1,dime(2)
      do i=1,dime(1)
           label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)

           Lx = i-xs(1)
           Ly = j-xs(2)
           Lz = k-xs(3)
           Lx = Lx - dime(1) * anint( Lx / dime(1) )
           Ly = Ly - dime(2) * anint( Ly / dime(2) )
           Lz = Lz - dime(3) * anint( Lz / dime(3) )
           

           checkx = (Lx)/real((R-.5d0))
           checkx = checkx*checkx
           checky = (Ly)/real((R-.5d0))
           checky = checky*checky
           checkz = (Lz)/real((R-.5d0))
           checkz = checkz*checkz
           check1 = checkx+checky+checkz
           

           checkx = (Lx)/real((R+.5d0));
           checkx = checkx*checkx;
           checky = (Ly)/real((R+.5d0));
           checky = checky*checky;
           checkz = (Lz)/real((R+.5d0));
           checkz = checkz*checkz;
           check2 = checkx+checky+checkz;
           
           if (check1.ge.1.d0.and.check2.lt.1.d0 ) then
              boundary_NP(i,j,k) = .true.
              cont_b = cont_b + 1
              lb2index(label) = cont_b
              temp_po(1,cont_b,ii) = i
              temp_po(2,cont_b,ii) = j
              temp_po(3,cont_b,ii) = k
           elseif (check1< 1.d0) then
              NP(i,j,k) = .true.
              cont_np = cont_np + 1
           endif
           
      enddo
   enddo
enddo

elseif (droplet) then

do k=1,dime(3)
   do j=1,dime(2)
      do i=1,dime(1)
         if(boundary(i,j,k).or.drop(i,j,k)) then
           label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)

           Lx = i-xs(1)
           Ly = j-xs(2)
           Lz = k-xs(3)

           checkx = (Lx)/real((R-.5d0))
           checkx = checkx*checkx
           checky = (Ly)/real((R-.5d0))
           checky = checky*checky
           checkz = (Lz)/real((R-.5d0))
           checkz = checkz*checkz
           check1 = checkx+checky+checkz
           

           checkx = (Lx)/real((R+.5d0));
           checkx = checkx*checkx;
           checky = (Ly)/real((R+.5d0));
           checky = checky*checky;
           checkz = (Lz)/real((R+.5d0));
           checkz = checkz*checkz;
           check2 = checkx+checky+checkz;
           
           if (check1.ge.1.d0.and.check2.lt.1.d0 ) then
              boundary_NP(i,j,k) = .true.
              if (boundary(i,j,k)) bnodes = bnodes -1
              boundary(i,j,k)=.false.
              drop(i,j,k)=.false.
              cont_b = cont_b + 1
              temp_po(1,cont_b,ii) = i
              temp_po(2,cont_b,ii) = j
              temp_po(3,cont_b,ii) = k
           elseif (check1< 1.d0) then
              NP(i,j,k) = .true.
              if (boundary(i,j,k)) bnodes = bnodes -1
              boundary(i,j,k)=.false.
              drop(i,j,k)=.false.
              cont_np = cont_np + 1
           endif
        endif   
      enddo
   enddo
enddo

endif 

NPnodes(ii) = cont_np 
NPbnodes(ii) = cont_b
enddo
End subroutine grid_NP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine find_normal1()
use type
use var 
implicit none
integer :: i, j, k, ii, jj,indx, label
real(dp) :: theta, phi, r, x,y,z
if (allocated(normal1)) deallocate(normal1)
if (allocated(tangent1)) deallocate(tangent1)
allocate(normal1(3,sum(NPbnodes)),tangent1(3,sum(NPbnodes)))
indx = 1
DO ii=1,number_NP
   Do jj=1, NPbnodes(ii)
      i = temp_po(1,jj,ii); j = temp_po(2,jj,ii); k = temp_po(3,jj,ii)
      label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
      lb2index(label) = indx
      x = real(i - NP_position(1,ii)); y = real(j - NP_position(2,ii)); z = real(k - NP_position(3,ii))
      ! appling periodic boundary
      if (x>dime(1)/2) then
         x=-dime(1)+x
      elseif (x<-dime(1)/2) then
         x=dime(1)+x
      endif
      if (y>dime(2)/2) then
         y=-dime(2)+y
      elseif (y<-dime(2)/2) then
         y=dime(2)+y
      endif
      if (z>dime(3)/2) then
         z=-dime(3)+z
      elseif (z<-dime(3)/2) then
         z=dime(3)+z
      endif

      r = sqrt(x*x + y*y + z*z)
      theta = acos(z/r)
      if (x/=0) then
         phi = atan(y/x)
      else
         phi = Pi /2.d0
      endif
      if (x<0.d0.and.y<0.d0) then
         phi = phi + Pi
      elseif (x<0.d0.and.y>0.d0) then 
         phi = phi + Pi
      elseif (x>0.d0.and.y<0.d0) then
         phi = phi + 2*Pi
      endif
      if (phi==0.d0.and.x<0) then
         theta = Pi - theta
      endif
      if (phi==Pi/2.d0.and.y<0) then
         theta = Pi - theta
      endif
      normal1(1,indx) = sin(theta)*cos(phi)
      normal1(2,indx) = sin(theta)*sin(phi)
      normal1(3,indx) = cos(theta)
      tangent1(1,indx) = cos(theta)*cos(phi)
      tangent1(2,indx) = cos(theta)*sin(phi)
      tangent1(3,indx) = -sin(theta)
      if (phi==0.d0.and.x<0) then
         normal1(:,indx) = -normal1(:,indx)
      endif
      if (phi==Pi/2.d0.and.y<0) then
         normal1(:,indx) = -normal1(:,indx)
      endif
      indx = indx + 1
   enddo
enddo
endsubroutine find_normal1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Subroutine anchoring_q()
use type
use var
implicit none
integer :: i, j, k, x, ii, jj, label, indx
real(dp), dimension(3) :: dr, normal, eigenvalues, r
real(dp) :: s
if (allocated(qb1)) deallocate(qb1)
allocate(qb1(6,sum(NPbnodes)))
indx = 1
Do ii=1,number_NP
   Do jj=1, NPbnodes(ii)
      
      if (NPtype=="xx") normal = tangent1(:,indx)
      if (NPtype=="zz") normal = normal1(:,indx)
    
      call director_tensor( Sbulk, normal, qb1(1:5,indx) )
      call eigen( qb1(:,indx), eigenvalues, r )
   
      s = maxval( eigenvalues )*3.d0/2.d0
      if ( abs(s - Sbulk) > 1.d-10 ) then 
      print *, '--> Problem with anchoring 1', i, r
      stop
      end if
      s = sqrt( dot_product( r,r ) )
      if ( abs(s - 1.d0) > 1.d-10 ) then 
      print *, '--> Problem with anchoring 2', i, r
      stop
      end if
      
      qb1(6,indx) = - qb1(1,indx) - qb1(4,indx)
      i = temp_po(1,jj,ii); j = temp_po(2,jj,ii); k = temp_po(3,jj,ii)
      label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
      q(:,label) = qb1(:,indx)
      indx = indx+1
   enddo
enddo
End Subroutine anchoring_q 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine find_boundary_neighbors1()
use var
use type
implicit none
integer :: i,j,k, indx, label
integer :: deltax,deltay,deltaz
if (allocated(boundary_neighbors1)) deallocate(boundary_neighbors1)
allocate(boundary_neighbors1(6,sum(NPbnodes)))
boundary_neighbors1 = 0
indx = 1

do k=1,dime(3)
   do j=1,dime(2)
      do i=1,dime(1)
         label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2) 
         if (boundary_NP(i,j,k)==.true.) then
            indx = lb2index(label)
            if (normal1(1,indx)>=0) deltax=1
            if (normal1(1,indx)< 0) deltax=-1
            if (normal1(2,indx)>=0) deltay=1
            if (normal1(2,indx)< 0) deltay=-1
            if (normal1(3,indx)>=0) deltaz=1
            if (normal1(3,indx)< 0) deltaz=-1
            boundary_neighbors1(1,indx) = 1 + (i-1+deltax) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
            boundary_neighbors1(2,indx) = 1 + (i-1+2*deltax) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
            boundary_neighbors1(3,indx) = 1 + (i-1) + (j-1+deltay)*dime(1) + (k-1)*dime(1)*dime(2)
            boundary_neighbors1(4,indx) = 1 + (i-1) + (j-1+2*deltay)*dime(1) + (k-1)*dime(1)*dime(2)
            boundary_neighbors1(5,indx) = 1 + (i-1) + (j-1)*dime(1) + (k-1+deltaz)*dime(1)*dime(2)
            boundary_neighbors1(6,indx) = 1 + (i-1) + (j-1)*dime(1) + (k-1+2*deltaz)*dime(1)*dime(2)
         endif
     enddo
   enddo
enddo
end subroutine find_boundary_neighbors1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1

End Subroutine make_np_channel

